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ABSTRACT 


An  examination  and  a  comparison  of  the  relative  merits  of 
the  finite  element  and  Fourier  series  methods  of  solving  radially 
loaded  circular  ring  problems  are  made.     The  procedure  employed 
to  evaluate  the  two  methods  is  to  use  each  method  to  solve  for  three 
different  load  conditions  and  to  compare  the  performance  of  the  two 
methods  on  the  basis  of  accuracy,    ease  of  usage,    and  equipment 
required.     The  results  indicate  a  satisfactory  accuracy  for  both 
methods  under  most  conditions.     The  Fourier  series  method  is 
superior  for  solving  problems  with  a  distributed  load  condition. 
The  finite  element  method  is  superior  for  solving  problems  with 
concentrated  loads. 
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CHAPTER  I 
INTRODUCTION 


The  structural  analysis  of  many  aerospace  vehicles,    such  as 
rockets,    re-entry  bodies ,   aircraft,    space  capsules,    etc.,    is 
frequently  accomplished  using  either  the  Fourier  series  method 
or  the  finite  element  method  in  conjunction  with  a  digital  com- 
puter.    A   considerable  amount  of  effort  has  been  devoted  to  the 
development  and  refinement  of  these  two  methods,    but  little 
attention  has  been  given  to  a  direct  comparison  of  the  methods. 
Fundamental  questions   regarding  the  similarities,    differences, 
relative  accuracy,    ease  of  usage,    and  advantages   or  disadvantages 
of  each  of  the  methods  have  not  been  answered  as   yet.      The  ob- 
jective of  this   study  was  to  seek  answers  to  some  of  these  funda- 
mental questions. 

The  procedure  employed  to  accomplish  the  objective  was  to 
apply  the  two  methods   of  analysis  to  a  portion  of  a  typical  aero- 
space structure  and  to  investigate  the  relative  merits  of  the  two 
methods.      The  portion  selected  for  analysis  was  a  thin  circular 
ring  with  a  rectangular  cross  section.      The  ring  was  selected 
because  it  represents,    in  a  simplified  manner,    the  circular 
cylinders  that  are  frequently  found  in  aerospace  vehicles.      Three 
different  load  conditions  were  applied  to  the  ring  for  each  of  the 
two  methods  of  analysis  so  that  the  effects  of  the  load  distri- 
butions on  the  performance  of  the  two  methods  could  be  investi- 
gated.     The  extent  to  which  the  analyses  were  carried  out  in 
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terms  of  numbers  of  significant  figures  carried,    numbers  of 
terms  evaluated  in  the  Fourier  series  method,    number  of 
elements  used,    and  use  of  double  precision  techniques  in  the 
FORTRAN  computer  programs,    was  greater  than  that  normally 
used  in  common  engineering  practice.     The  reason  for  such 
depth  was  to  bring  out  any  subtleties  peculiar  to  either  of  the 
two  methods  that  might  not  appear  if  normal  engineering    accu- 
racy were  used. 

The  author  gratefully  acknowledges  the  guidance  and 
assistance  of  Professor  Robert  E.    Ball  in  the  preparation 
of  this  thesis. 
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CHAPTER  II 

DERIVATION  OF  THE  GOVERNING  EQUATIONS 
FOR   THE  FOURIER  SERIES  ANA  LYSIS 


Consider  the  thin  ring  shown  in  Figure   1  where  a  is  the  radius 
of  the  ring  middle  surface;   b,    the  ring  width;    and  t,    the  thickness. 
The  angle  (J)  is  the  generator  of  the  ring  and  is  positive  in  a  counter 
clockwise  direction.     All  points  in  the  ring  are  located  by  the 
coordinates   (J)  and  z,    where  z  originates  at  the  middle  surface  and 
is  positive  in  the  outward  radial  direction.      The  normal  force  N^> 


FIGURE   1 
RING  CONFIGURATION 
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the  shear  force  Q^,    the  bending  moment  M.  ,    and  the  applied 
radial  loading  Pr  are  shown  in  the  positive  direction.      The  dis- 
placements,   also  shown  in  the  positive  direction  in  Figure    1,    are 
w  in  the  radial  direction  and  v  in  the  tangential  direction. 

Assuming  that  the  ring  retains  its  circular  shape  after  de- 
formation,   the  equations  of  equilibrium  can  be  written  as 

aQ*   -      =     0  (2-la) 

a  4> 


and 


<4  4> 


^ 


(2-lb) 


"    Q^    =       0  (2-lc) 


Using  Equation  (2-la)  to  eliminate  Q^  from  Equations  (2-lb)   and 
(2-lc)  leads  to  the  two  equilibrium   equations 

i$  a  4)  (2-2a) 


2D  ^  d> 


d  (J)z  (2-2b) 

The  circumferential  strain,  £  ^     ,     at  any  distance  z  from 
the  middle  surface,    is   given  in  Eeference   1  as 

c    -     ° £_       +    (2-3) 

•        adcf)         a.  fa  +  2.)    d^*  a  +  z. 

This   relationship  is  based  upon  the  assumption  that  points  on  a 
normal  to  the  middle  surface  prior  to  deformation  remain  on  the 


16 


normal  after  deformation  and  that  the  thickness  of  the  ring 
remains  constant.     Substituting  Equation  (2-3)  into  the  plane 
stress  form  of  the  elastic  constitutive  law  leads  to 


<r«>  =  E 


1 


\r 


1        cTuj 


4- 


UU 


a  +  z. 


I  a<i$         a.    (a  +  z-)     \<\>i 

where  (Ti  is  the  circumferential  stress  at  the  location  (({),    z) 
and  E  is  the  modulus  of  elasticity. 

The  normal  force  and  bending  moment  are  given  by 


(2-4) 


N«,    = 


I 


J      1 


(T^dz. 


(2-5a) 


and 


M*=jt*  (r*zd 


X 


(2-5b) 


Substituting  Equation  (2-4)  into  Equations  (2-5a)and  (2-5b), 
integrating  and  applying  the  limits  yield 


M  & 


d$ 


-t-  u; 


+ 


K 


A1^ 
d<&1 


-r-  UU 


(2-6a) 


and 


M$~    a* 


+  uu 


(2-6b) 


w 


here  D  =  Et  and  K  =  Et3/12  .        * 


*In  order  to  use  this  form  of  D  and  K  the  natural  logarithmic 
terms  that  result  from  the  integration  of  Equations  (2-5a)  and 
(2-5b)  must  be  expanded  in  a  series  in  terms  of  t/2a  where  terms 
to  the  fifth  power  and  above  are  dropped. 
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The  problem  has  now  been  reduced  to  one  with  four  utikm  wns, 
v,    w,    N,k,   and  M^,    and  four  equations;     two  of  the  equations  are 
obtained  from   the  conditions  of  equilibrium 


a 


d<t> 


d$ 


o 


(2-2a) 


aM^-a2^  + 


4> 


=  o 


(2-2b) 


and  two  are  obtained  from  the  application  of  the  elastic  consti- 
tutive law  and  the  strain-displacement  relationship 


do 


4-  uu 


K 


^3    L  <i  4>^ 


uu 


(2-6a) 


"V* 


azo) 


+  uu 


(2-6b) 


I    d<^ 

Taking  the  appropriate  derivatives  with  respect  to  (J)  of  Equations 
(2-6a)  and  (2-6b)  and  substituting  the  results  into  Equations  (Z-2a) 
and  (Z-2b)  yield 

d  r  Av 


d(()  L    d<^ 


-+-  uo 


and 


dvr 
d$> 


W  +A 


=  O 


d  u;        ~d  tu 


(2-7a) 


La  <j>h- 


d(p: 


(JL) 


D 


=  O 


(2-7b) 


2  /         2  ' 
where  k  =  t  /12a    .     Equations  (2-7a)  and  (2-7b)  are  the  governinj 

differential  equations  for  a  thin  ring  under  the  assumptions  that 
have  been  made. 
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CHAPTER  in 

DEVELOPMENT  OF  THE  FOURIER  SERIES  SOLUTION 
TO  THE  GOVERNING  EQUATIONS 


The  starting  point  in  the  development  of  the  Fourier  series 
solution  to  Equations  (2-7a)  and  (2-7b)  is  to  assume  that  the 
applied  load,    Pr  ,  is  a  general  function  of  (J>,    symmetrically  dis- 
tributed about  two  perpendicular  planes,    one  through  (j)  =   0  and 
0  -  1"f\   and  one  through  0=-'TT/2.     Thus,    the  load  can  be  repre- 
sented by  the  Fourier  cosine  series      * 
oo      „ 

?X=L    r w  c°s  V"  <t>  (3-1) 

where  the  coefficients  Prn  are  given  by 

P0=  t/    P,  H  (3"2a) 

^o 
and 

^=    ~TM       f>rCoSm*c^<f>  (3_2b) 

The  radial  displacement,    w,    which  is  an  even  function  of  d) 

for  the  ring  under  the  load  distribution  considered,    is  also  written 

in  a  Fourier  cosine  series  as 

oo 
W  =E     UUmCOS    W\    (J)  (3-3) 


*  The  symmetry  conditions  on  Pr  lead  to  the  elimination 
of  all  odd  values  of  m,     i.  e.  ,    Pm  =  0  for  all  odd  values  of  m. 
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The  tangential  displacement  v  is  an  odd  function  of  (J)  and  is 
given  by 


oo 


\f  =  I    ^m     Mn    W(j)  (3-4) 

The  governing  differential  equations  are  used  to  obtain  w 
and  v       as  a  function  of  Pm.     For  the  axisymmetric  mode, 
where  m  =  0,    the  substitution  of  Equations  (3-1)  and  (3-3)  into 
Equation  (2-7b)  yields 


UU0  =     — ^—  (3-5) 


For  the  general  case  of  m  >0,    substitution  of  Equations  (3-3) 
and  (3-4)  into  Equation  (2-7a)  results  in 

mx  vm  +  m  u)m  =  o  (3-6) 

Substituting  Equations  (3-1),    (3-3),    and  (3-6)  into  Equation  (2-7b) 

gives 

^  f 


ID      - 


Therefore,  according  to  Equation  (3-6) 

Substituting  these  expressions  for  w       and  vrn  into  Equations   (3-3) 
and  (3-4)  leads  to. 

U,g   -gl^-   +  -«lf      Pwcosw>»  (3.?a) 
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v  =  - y    Laa (3-7b> 

All  other  quantities  of  interest  are  determined  by  using   these 
expressions  for  w  and  v  in  the  appropriate  equations.     For 
example,    according  to  Equations  (2-6a)  and  (2-6b),    the  normal 
force  is  given  by 


p  ~  ?m    COSM 


(3-8) 


And  the  bending  moment  is 


W*=-TT^--aE  ^ (3-9> 

To  permit  a  direct  comparison  between  the  Fourier 

series  and  finite  element  analyses  of  the  ring,    w,    Ni,   and  Mi, 

are  selected  for  evaluation  by  both  methods  at  the  points  (j>  ~  0, 

30,    60,    and  90  degrees.     In  Chapter  IV,    P      ,    and  the  specific 

m 

equations  for  w,    N ,  ,   and  M^  are  developed  and  evaluated  at 
these  values  of  ty  for  three  different  load  conditions. 
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CHAPTEE   IV 
APPLICATION  AND  EESULTS  OF  THE  FOUEIEE  SEEIES  SOLUTION 

The  three  loads  selected  for  use  in  this  study  are  shown  in 
Figure  2.     All  three  loads  are  radial  loads,    directed  inward,    and 
are  symmetrically  distributed  about  two  planes,    one  through  0=0 
and  (1)   =  TT    and  one  through  0  =  tTT/2.      Case  1  is  a   "line"  or  "knife 
edge"  load  of  P  pounds  distributed  over  the  width  of  the  ring  at  0  =  0 
and  0  =  TT- 


P(lb 


V(\V$\ 


Case   1             oc  z 

0  degrees 

Case  2            o<   z 

30  degrees 

Case  3            *^  = 

60  degrees 

FIGUEE  2 

LOADS  CONSIDEEED 
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Case  2  and  Case  3  loads  are  uniform  pressure  loads,    each  with  a 
total  value  of  P  pounds  covering  the  areas  shown  in  Figure  2.     The 
uniform  loads  are  initially  represented  as  P  pounds  distributed 
over  the  angle  2°<  so  that  the  surface  load  P     is  given  by 

P 


P  =-   - 

r  Zc*  4  b 

Thus,    substituting  Equation  (4-1)  into  Equation  (3-1)  gives 


(4-1) 


P 


oo 


Z<*ab 


-    -I       L  C03   **  * 


m  =  0,1,4 


Accordingly  Equation  (3-2a)  becomes 


l  =  - 


Z<*a  b 


so  that 


K~- 


(4-2) 


TTab 

Similiarly,    substituting  Equation  (4-1)  into  Equation  (3-2b),    inte- 
grating and  applying  the  limits  yield 

P 


m  f  a  b  m<=< 


[sin  m=<  -  stm  mft— °<)] 


Therefore 


z  P  s\n  m°^ 

™  TT  a  b   vn  o< 

for  m  even.  For  odd  values  of  m 


(4-3) 


as  originally  noted  in  Equation  (3-1). 
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Applying  Equation  (4-2)  and  (4-3)  to  Equations  (3-7a),    (3-8) 

and  (3-9)  leads  to  the  Fourier  series  solution  for  w,    N     and  M^ 

tf>  <p 


uu=  - 


a? 


00 


2iaP      ^      Mn    W\  o<.  c  o  S  W\  $ 


(4-4a) 


N* 


+   V    

ttVd       trb  °<  m  ( mx-i) 


(4-4b) 


and 


d-tlP  £qP     y       MY\  m«x     to±VY\§ 


U*       ma+Jk)     irb 


(4-4c) 


The  non-dimensional  forms  of  these  three  equations,    denoted  by 
w,    Nx  and  M,k,   and  obtained  by  using  Jk  —  t    /12a    ,   D=  Et,    and 


> 


^ 


I  =  bt    /12,    are 


(jj 


ul 


U) 


a>p 


ir^Ad-e^ 


oo 

■I 


sm  w\°<  c  os  w\(t> 


(4-5a) 


NL  = 


kN* 


1  S2,        S  I  Y\    VY\<*-     C  OS  W  (J> 


(4-5b) 


and 


M. 


ira^  K(ui)       ir 


S  \  Y\  W  <*•    COS    VV^ 


(4-5c) 
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These  expressions  are  applicable  only  to  Case  2  and  3  loads. 
For  Case  1  load,    c<  — ►»  0.     Thus,    taking  the  limit  of  Equations 
(4-5a),    (4-5b)  and  (4-5c)  as   o<  — >- 0    leads  to  the  solution  for 
Case   1  load 


X                          %     ?2               COS  W\$ 
Co  = 1        ; rr  (4-6a) 


N £ -i-x    --f-  (4-6b) 

and  TY\=i,^ 


M     -_    . y     (4-6c) 

m  =  2.,  t 

The  term   -  I/TT  a       A   (1+  k)  that  appears  in  w  and  Mx  will  now 
be  examined  for  its  contribution  and  significance.      The  expressions 


iU^-jbM       fl«i      <}  IV 


M*a(<M 


give  the  axial  strain  energy  and  the  bending  strain  energy  respec- 
tively of  an  element  bt<x(d(|))  of  the  ring.  They  show  that  the  axial 
strain  energy  is  proportional  to   1/A  and  the  bending  strain  energy 

is  proportional  to  1/1.      Using  U  ~   1/A  and  U,~    1/1  leads  to 

a  b 

U    /Ul~  I/A.     In  Reference  2,    the  axial  strain  energy  for  such  a 
curved  element  is  shown  to  be  very  small  compared  to  the 
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bending  strain  energy  of  the  element  provided  that  t/a  is  small.  * 
Thus,    the  term   -  I/TT  a       A   (1  +  k)  in  w  and  M,  is  related  to  the 
extension  of  the  middle  surface  of  the  ring  and  should  be  small. 
If  the  middle  surface  is  considered  inextensional,    the  axial 
strain  energy  of  an  element  of  the   ring  is   zero,  and  the  term 
-I/TT    a       A(l+-k),    which  is  proportional  to  U    /Uv  ,    does  not 
appear  in  the  expressions  for  w  and  M..     In  order  to  determine 
the  magnitude  of  the  contribution  of  the  axial  strain  energy  to  the 
solution,   a  t/a  ratio  must  be  selected.     As  an  example,   a  ratio 
of  t/a  =  1/10  is  chosen  and  w  and  M.   are  calculated  with  and 
without  the  term  containing  I/A  at  (I)  =  0,    30,    60  and  90  degrees 
for    all  three  load  conditions.     The  results  of  these  calculations 
are  found  in  Table  I,    where  the  values  referred  to  as  "EXT" 
include  the  -I/TT  a2  A   (1  +  k)  term  and  the   "NON-EXT"  values   do 
not.     A  s  is  evident  from   Table  I  the  differences  between  the  "EXT" 
and  "NON-EXT"  values  are  very  small,    or  nonexistent  for  the 
accuracy  retained,    when  t/a  a  10. 

An  important  consideration  in  the  application  of  the  Fourier 
series  method  is  the  number  of  terms  of  the  series  required 
to  give  a  desired  accuracy  in  the  solution,     i.  e.  ,    how  fast  does 


#This  statement  assumes  there  are  no  distributed  loads  on 
the  element,    i.e.,    only  concentrated  loads  are  considered.     This 
thesis  considers  distributed  loads,   and  thus  the  statement  is  not 
completely  applicable. 
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TABLE  I 
FOURIER  SERIES  EXTENSIONS  L  AND  NON  -  EXT  ENSIONAL  RESULTS 


CASE  1 

LOAD 

CASE 

2  LOAD 

CASE  3 

LOAD 

0 

Ext. 

Non-ext. 

Ext. 

Non-ext. 

Ext. 

Non-ext, 

W 

-0. 0747 

-0.O744 

-0. 0599 

-0.  0596 

-0. 0289 

-0.  0287 

0° 

fi, 

0.  0 

0.  0 

-0.  128 

-0.  128 

-0.239 

-.0.  239 

M$ 

0.  318 

0.  318 

0.  190 

0.  190 

0.  079 

0.  080 
-0.  0149 

w 

-0. 0337 

-0.  0334 

-0. 0289 

-0. 0287 

-0.  0152 

30° 

s» 

-0.  250 

-0.  250 

-0.  239 

-0.  239 

-0.271 

-0.  271 

M* 

0.  068 

0.  068 

0.  079 

0.  080 

■0.  047 

0.  048 

w 

0.  0361 

0.  0364 

0.  0295 

0.  0298 

0. 0141 

0.  0143 

60° 

s» 

-0.  433 

-0.433 

-0.413 

-0.  413 

-0.  358 

-0.  358 

M+ 

-0.  115 

-0.  115 

-0.  095 

-0.  095 

-0.  040 

-0.  040 
0.  0298 

w 

0.  0680 

0.  0683 

0.  0571 

0.  0574 

0.  0295 

90° 

1> 

-0.  500 

-0.  500 

-0.477 

-0.  477 

-0.413 

-0.413 

_ 

0 

-0.  182 

-0.  182 

-0.  159 

-0.  159 

-0.  095 

-0.  095 
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the  series  converge?     This  information  for  the  inextensionai 
form  of  the  solution  is   shown  in  Tables  II,    III  and  IV  for  the  first 
4  terms  of  w,    N^and  M-  at  0  =  0,    30,    60  and  90  degrees  for  all 
three  load  conditions.    The  number  of  terms  required  for  con- 
vergence of  a  series  to  four  decimal  place  accuracy  for  the 
displacement  and  three  decimal  place  accuracy  for  the  forces 
and  bending  moments  is  shown  in  the  last  column  with  the  con- 
verged value  of  the  series. 
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TABLE  II 


FOURIER  SERIES  NON 

FOR 

-EXTENSIONAL  CONVERGENCE 
CASE  2  LOAD 

Term 

i 

2 

3 

4 

Converged 
Value 

i 

♦ 

m 

2 

4 

6 

8 

&m 

-0. 0707 

-0.  0736 

-0.  0741 

-0. 0742 

-0.  0744  at  m  = 

14      j 

0° 

S<m 

-0.  106 

-0.  064 

-0. 045 

-0. 035 

0.  000   at  m    ==. 

534 

&*m 

0.  212 

0.  255 

0.  273 

0.  283 

0.  318     at  m  = 

1 

394 

<—  m 

-0.  0354 

-0.  0340 

-0. 0334 

-0.  0334  at  m  = 

6 

30° 

*—  m 

-0.  212 

-0.  233 

-0.  252 

-0.257 

-0.  250     at  m  = 

34 

^  m 

0.  106 

0.  085 

0.  067 

0.  062 

0.  068     at  m   = 

52 

(—>  m 

0.  0354 

0.  0368 

0. 0363 

0.  0363 

0. 0364  at  m   = 

14 

60° 

Sm 

-0.  424 

-0.  446 

-0.  427 

-0.432 

-0.  433     at  m   = 

24 

S*m 

-0.  106 

-0.  127 

-0.  109 

-0. 114 

-0.  115     at  m   = 

38 

Sm 

0.  0707 

0. 0679 

0. 0684 

0. 0683 

0.  0683  atm    = 

8 

90° 

s 

0  m 

-0.  531 

-0.  488 

-0.  506 

-0.496 

-0.  500     atm    r 

26 

a  m 

-0.  212 

-0.  170 

-0.  188 

-0.  178 

-0.  182     atm    = 

38 
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TABLE  III 

FOUEIEE  SEEIES  NON-EXTENSIONA  L  CONVEEGENCE 

FOE  CASE  2  LOAD 


Term 

1 

2 

3 

4 

Converged 
Value 

♦ 

m 

2 

4 

6 

8 

0° 

D"m 

-0.  058 

-0.  143 

0.  175 

-0. 0547 
-0.  125 
0.  193 

-0.  0597 
-0.  125 
0.  193 

-0. 0596 
-0.  127 
0.  191 

-0.  0596  at  m  =     8 

-0.  128     at  m  =   10 

0.  190     at  m  =  20 

o 
30 

-0. 0292 
-0.  231 
0.  088 

-0.  0287 
-0.239 
0.  079 

-0.  239 
0.  079 

-0.  238 
0.  080 

-0.  0287  at  m  =     4 

-0.  239     at  m  =  10 

0.  080     at  m  =  14 

60° 

0.  0292 
-0.  406 
-0.  088 

0. 0298 
-0.415 
-0.  097 

-0.415 
-0.  097 

-0.414 
-0.095 

0.  0298  at  m  =     4J 
-0.  413     at  m  =  44 
-0.095     at  m  =     8 

90° 

£™m 

0.  0585 
-0.  494 
-0.  175 

0.  0573 
-0.476 
-0.  158 

0.  0573 
-0.  476 
-0.  158 

0.  0574 
-0.478 
-0.  160 

0.  0574  at  m  =     8 
-0.  477     at  m  =  10 
-0.  159     at  m  =  10 

■■■■'—■■■■                  '                                       ■■■»■■■■■■■■ 
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TABLE  IV 

FOURIER  SERIES  NON -EXTENSIONS  L  CONVERGENCE 

FOR  CASE  3  LOAD 

Converged 

Term 

1 

2 

3 

4 

Value              i 

! 

0    Itn                   2 

4 

6 

8 

i 

i 

i 

! 

i 

-     m 

-0.  0292 

-0. 0287 

-0.  0287  at  m  =    4; 

0° 

?-"m 

-0.  231 

-0.239 

-0.  239 

-0.238 

-0.  239     at  m  =  10; 

^Mm 

0   088 

0.  079 

0.  079 

0.  080 

0.  080     at  m  =  14! 

m 

-0. 0146 

-0. 0149 

-0.  0149  at  m  =    4 

30° 

"Nm 

-0.  274 

-0.  270 

-0.  270 

-0.  271 

-0.  271     at  m  =     8 
0.  048     at  m  =  14 

i 

rr> 

0.  044 

0.  048 

0.  048 

0.  048 

I^m 

0. 0146 

0.  0143 

0.  0143  at  m  =    4 

60° 

<^Nro 

-0.  362 

-0.  358 

• 
-0.  358     at  m  =    4; 

7^m 

-0. 044 

-0. 039 

-0. 039 

-0.  040 

i 
-0.  040     at  m  =     8 

■ 

0.  0292 

0.  0298 

0.  0298  at  m  =    4 

»0° 

DS« 

-0.  406 

-0.415 

-0.415 

-0.414 

-0.  413     at  m  =  44 

J 

-0. 088 

-0.  097 

-0.  097 

-0.095 

-0.  095     at  m  =     8 
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CHAPTEE  V 
DEVELOPMENT  OF  THE  FINITE  ELEMENT  SOLUTION 

The  finite  element  method  is  based  upon  the  concept  of 
replacing  the  actual  continuous  structure  with  a  number  of  struc- 
tural elements  of  finite  size,    leading  to  an  assembled  structure. 
The  requirements  of  equilibrium  and  continuity  at  the  element 
joints  or  nodes  lead  to  a  set  of  simultaneous  algebraic  equations. 
There  are  two  general  methods  for  formulating  these  equations, 
the  displacement  method  and  the  force  method.     The  displace- 
ment method,    also  known  as  the  stiffness  or  direct  stiffness 
method,    is  the  one  used  in  this  thesis  to  solve  the  ring  problem. 
This  method  is  discussed  in  detail  in  Reference  3.     A  brief  ex- 
planation of  the  method  is  given  here  for  continuity. 

In  the  displacement  method  the  displacements  and  rotations 
at  the  nodes  of  the  elements  comprising  the  structure  are  taken 
as  the  unknown  quantities.      The  matrix  equation  which  relates 
the  external  forces  to  the  displacements  of  an  assembled  structure 
is 

\x]    -     K      [n\     4-  (x°l  (5-1) 

Matrix   *.  X^    is  the  column  matrix  of  nodal  forces,    i.  e.  ,    forces 
that  act  upon  the  nodes  of  the  element.     In  a  general  three  dimen- 
sional structure  these  may  include  bending  and  twisting  moments 
as  well  as  rectilinear  forces.     The  matrix     K      is  a  square 
symmetric  matrix  for  the  entire  structure.     This  matrix  is 
known  as  the  stiffness  matrix  and  is  assembled  from  the  individual 
stiffness  matrices  of  the  elements  forming  the  structure.      The 
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entries  in  the  element  stiffness  matrices  are  called  stiffness 
influence  coefficients.     They  are  determined  by  displacing  a 
node  of  the  element  by  one  unit  in  the  direction  of  one  of  the 
degrees  of  freedom  permitted  for  that  node  and  calculating 
the  nodal  forces   required  to  maintain  that  unit  displacement  and 
to  prevent  movement  in  any  other  degree  of  freedom  for  the 
nodes  on  the  element.      The  amount  of  force  required  at  a  node 
to  maintain  this  displacement  state  is  one  stiffness  influence 
coefficient.     The  complete  set  of  stiffness  influence  coefficients 
is   determined  by  repeating  this  procedure  for  each  degree  of 
freedom  at  every  node.      The  column  matrix  ^u]  is  the  matrix  of 
displacements  and  rotations  which  may  take  place  at  a  node.      The 
column  matrix  JX    j  is  a  matrix  of  initial  forces.      They  are  deter- 
mined by  clamping  all  of  the  nodes   of  the  assembled  structure 
and  calculating  the  external  force  required  at  each  node  to 
maintain  this  zero  displacement  state  under  the  applied  load  con- 
dition.     Thus,    loads  distributed  between  nodes  and  concentrated 
loads  not  located  at  nodes  can  be  treated. 

There  are  various  techniques  available  for  the  manipulation 
and  solution  of  Equation  (5-1).     The  one  used  in  this  study  is  to 
transpose  (x°)  to  the  left-hand  side  of  Equation  (5-1).     This  leads 
to  a  single  column  matrix  of  forces   |{X|  -       (X    \\     .     Next, 

the  matrix  equation  is   reduced  in  order  by  a  consideration  of  the 
boundary  conditions  on  the  structure.     If  the  boundary  conditions 
specify  zero  displacement  at  some  nodes,    then  those  nodal  dis- 
placements are  no  longer  unknowns  and  are  removed  from     |u  \ 


along  with  the  row  of  the  equation  and  the  column  of 
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K 


that 


correspond  to  that  zero  displacement.     The   remaining  matrix 
equation,   which  is  referred  to  as  the  reduced  matrix  equation  and 
is  subscripted  with/",     is 


The  unknowns  in  Equation  (5-2)  are  contained  in  the  matrix  of 
displacements  <un|  .      Equation  (5-2)  is  a  set  of  simultaneous 
linear  algebraic  equations     that  is  solved  on  a  digital  computer 
using  the  FORTRAN  IV  subroutine  "DSIMQ,  "  a  description  of 
which  appears  in  Appendix  C.     The  outputs  of  "DSIMQ"  are  the 
nodal  displacements    \u-p\    . 

In  order  to  determine  the  nodal  forces    ^X\    ,    the  stiffness 
matrix  equation  is  solved  for  each  element.     These  equations, 
subscripted  in   ,    are 

!x4  =  [K<\]Ki+!xV)  (5_3) 

The  unknowns  in  this  equation  are  in    \  XyJ    since  the  displacements 
{  uA    are  obtained  from  the  appropriate  locations  in      )UP     and   |u|, 
The  matrix    JX^  is  obtained  by  adding  the  matrix     jX^.  j  ,     which 
is  determined  by  selecting  values  from  appropriate  locations  in 

\X        ,    to  the  matrix    iKu!     )uy,l       .       With    Wy^   and     JX^J  known, 
all  of  the  nodal  displacements  and  forces  have  been  determined, 
and  the  problem  is  solved  since  the  forces  and  displacements  at 
any  interior  location  of  an  element  can  be  determined  from  the 
applied  load  and  nodal  forces  and  displacements. 
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Applying  the  J  Lnite  element  method  to  the  ring  problem  under 
consideration  requires  the  selection  of  the  element  geometry. 
Either  straight  or  curved  elements  can  be  used.     Curved  elements 
are  selected  since  they  more  naturally  match  the  shape  of  the  ring 
and  result  in  displacements  and  forces  that  are  readily  interpreted 
without  requiring  a  change  in  coordinates.      The  curved  element 
used  is  drawn  in  Figure  3  with  the  nodal  forces  N,    Q,    N,    Q, 
bending  moments  M .  ,    M_  ,    displacements  w,    v ,  ,    w,    v?  and 
rotations   9.,    0_  shown  in  the  positive  sense.      The  element  covers 
the  angle    fi     .       A   stiffness  matrix  for  this  element  has  been 
derived  in  Eeference  3*     based  on  the  assumption  that  the  thick- 
ness   of  the  element,    t,    is  small  compared  to  the  radius  of 


"ij^l 


\ 


Q*>-*U 


T*a*ei 


Qirwi 


FIGURE  3 
FINITE  ELEMENT  CONFIGURATION 


*A   more  detailed  derivation  and  a  comprehensive  study  of 
the  curved  element  appear     in  Reference  4. 
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curvature,    a,    so  that  only  the  strain  energy  of  bending  need  be 
considered  during  the  derivation.      This  assumption  neglects  the 
shear  distortion  energy,    the  energy  resulting  from  any  coupling 
between  bending  moment  and  normal  force,    the  energy  of  axial 
deformation,    and  the  displacements  that  correspond  to  these 
three  forms  of  strain  energy.     The  stiffness  matrix  equation  from 
Eeference  3  for  the  curved  element  can  be  given  in  a  form  corre- 
sponding to  Equation  (5-2)  as 


Nl 
Q, 


Mj/a 


r. 


-w 


1 


Nl° 


aO 


{**.     =5  N 


Q 


1 


Q- 


V  Mz/a; 


■w2 

a© 


2/ 


M^/a 


Q 


V    (5-4) 


M2    /a  } 


where      K^   is  the  six  by  six  matrix  of  stiffness  influence 
coefficients  for  the  curved  element.     The  matrix      K^    ,    each 
term  of  which  is  a  function  only  of  (~>  ,     will  not  be  written  out  in 
functional  form  because  it  is  quite  lengthy.      For  this  thesis  a 
double  precision  digital  computer  subroutine,    "RELM,"  was 
written  which  computes      K«    for  an  input  of  [~>  .     A   copy  of 
"EELiM"  is  included  in  Appendix  C.     The  stiffness  matrix       K^ 
computed  by  "EELM"  for  /3    =    7.5,     15,   and  30  degrees  is 
presented  in  Tables  V,   VI  and  VII  respectively. 

With       Kf,    determined,     K    can  be  assembled.     Equation  (5-1) 

is  completed  for  the  ring  by  using  the  loading  conditions  to  obtain 
sX    |  and  the  appropriate  portions  of     <  X(    .     Because  of  the 
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TABLE  V 


K, 


for   ft   =  7.  5  d 


egrees 


18,665,685      1,223,063     26,675      -18,665,640       1,223,762     -26,721 


1,223,063  85,518        2,099         -1,223,762  74,855        -1,400 


26,675  2,099 


69 


26,721  1,400 


23 


^18,665,640    -1,223,762-26,721         18,665,685     -1,223,063        26,675 


1,223,762  74,855        1,400         -1,223,063  85,518        -2,099 


-26,721  -1,400 


■23 


26,675  -2,099 


69 


3,  306 


76, 035 


3,  329 


TABLE  VI 
K^  for    p  =  1 5    degrees 


524 


34 


3,  329 


-350 


11 


3,  306 


576,897     75,862   3,306     -576,874      76,035  .-3,329 


75,862     10,658     524      -76,035       9,339     -350 


340 


524 


11 


-576,874     -76,035  -3,329      576,897     -75,862    3,306 


9,339     350      -75,862      10,658     -524 


34 
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TABLE  VII 
[k^I         for    ft    =30  degrees 


17,243  4,598  399  -17,232                  4,640            -411 

4,598  1,317  130             -4,640 

399  130  17                  -411 

-17,323  -4,640  -411              17,243 

4, 640  1,  158  87             -4, 598 

-411  -87  -6                   399 


1,  158 

-87 

87 

-6 

4,  598 

399 

1,  317 

-130 

-130 

17 
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symmetrical  nature  of  the  applied  loads  only  one  quadrant  of  the 
ring  needs  to  be  considered.     The  quadrant  from   (J)  =  0    to  90  degrees 
is   selected  so  that  a  direct  comparison  with  the  Fourier  series 
solution  for  w,    N.   and  M,   can  be  made  for  0=  0,    30,      60  and 
90  degrees . 


' 
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Chapter  vi 
application  and  results  of  the  finite  element  solution 

In  order  to  apply  the  finite  element  solution  to  the  ring 
quadrant  (J)  —  0  to  90  degrees,    Equation  (5-1)  is  assembled  for  the 
entire  quadrant  using  Equation  (5-4)  for  each  of  the  elements  com- 
prising the  quadrant.     KJ    matrices  for    (->   —  7.  5,    15  and  30  de- 
grees are  used  in  this  thesis  so  that  the  quadrant  is  divided  into 
12,    6  and  3  elements   respectively  for  the  solution  of  each  of  the 
three  load  conditions.      This  assembled  equation  is   reduced  and 
solved  in  the  manner  described  in  Chapter  V.      The  process  of 
assembling,    reducing,    and  solving  Equation  (5-4)  is  best  described 
by  using  an  example.      Case  2  load  with  /->     =    30  degrees  is   selected 
as  the  example. 

Figure  4  shows  the  assembled  quadrant  and  the  node  numbering 
sequence  for  the  three  elements  in  the  quadrant.      For  three 


FIGURE  4 
3  ELEMENT  QUADRANT  FOR  CASE  2  LOAD 
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elements     |K        is  a    12    x  12  square  symmetric  matrix  which  is 
assembled  using  three  of  the  6x6       KfJ  matrices  from  Table  VII. 
When  assembled  for  this  example,    Equation  (5-1)  appears  as 


N 


Qj  =   0 
iV^/a 


N2   =   0 
Q2  =   0 


f  v'c=   0 


-w 


a0     =   I 


-w- 


M7/a  =  0      EI      K    12  x  12     J   a0o 
N3   =   0 


\ 

V  5 

Qi° 

M^/a 

o 
N2 

Q2° 

1 

M2°/a^ 

Q3  -0 

M3/a 

N„ 


Q4  =  o 

M4/a 


; 


•VI 


3 
a9o 


-w 


a94=0; 


N- 


Q3° 

M3°/a 

N4° 


Q< 


M4   /a 


(6-1) 


w 


here  the  zero  entries  in  the     |X\    and     iu  I     matrices  are 
obtained  from  the  symmetry  characteristics  of  the  load.      The 
initial  forces  in  \X   ]  are  derived  in  Appendix  B  and  presented  in 
Table  BI.      They  are  transposed  to  the  left-hand  side  of  Equation 
(6-1)  to  give 


41 


/  N 


M2/a 
0 
0 
0 


.  12709678^ 
. 22182251 
. 01944670 
. 12707678 
. 22182251 
.01944670 


i     0 
\ 
0 

0 

N 
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The  reduced  matrix  equation  is  created  from  Equation  (6-2)  by- 
striking  out  the  rows  of  the  equation  and  columns  of         K 
12  x  12.      The  reduced  equation  with        Kp     an  8x8  matrix 
appears  in  Table  VIII. 

A   computer  program,    "BETA    30,  "  a  copy  of  which  appears 
in  Appendix  C,    solves  Equation  (6-3)  for      \up\  by  using  the 
subroutine  "DSIMQ.  "     For  this  example  the  output  of  MDSIMQ"     is 
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TABLE  IX 
FINITE  ELEMENT  EESULTS  FOR   CASE  1  LOAD 


/3  =  7.  5° 


/3=  15° 


ft    =30° 


30 


60^ 


w 


N 


M 


w 


N 


M 


w 


N 


M. 


•0.  0744 
0.  000 
0.  318 


0.  0334 
0.  250 
0.  068 


0. 0364 
-0.433 
-0.  115 


•0.  0744 
0.  000 
0.  318 
■0.  0334 
■0.  250 
0.  068 


W 

0.  0683 

0. 0683 

-0.  0683 

90° 

N 

!  -0.  500 

• 

-0.  500 

-0.  500 

..  .       _J 

M 

* 

I  -0.  182 

-0.  182 

-0. 182 

0. 0364 
■0.  433 
■0.  115 


-0.  0744 

0.  000 

0.  318 

-0.  0334 

-0.  250 

0.  068 


0. 0364 
-0.433 
-0.  115 
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TABLE  X 
FINITE  ELEMENT  RESULTS  FOE  CASE  2  LOAD 


*. 

&    =    7.  5° 

/3  =  15° 

/3     =  30° 

w 

-0.  0596 

-0. 0596 

-0.  0596 

0° 

N 

-0.  128 

-0.  128 

-0. 128 

^ 

0.  190 

0.  190 

0.  190 

w 

-0. 0287 

-0. 0287 

-0.  0287 

30° 

*♦ 

-0.239 

-0. 239 

-0.  239 

■;  M* 

0.  080 

0.  080 

0.  080 

.  w 

0.  0298 

0. 0298 

0.  298 

60° 

Nv" 

-0. 413 
-0.  095 

-0.413 
-0.  095 

-0.  413 
-0.  095 

w 

0.  0574 

0.  573 

-0.  0573 

90° 

S* 

-0.  477 

-0.477 

-0.  477 

M^ 

-0.  159 

-0.  159 

-0.  159 
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TABLE  XI 
:  ELEMENT  RESULTS  FOR  CASE  3  LOAD 


♦ 

/3  =  7.5° 

#    =  15° 

/S  =  30°    ! 

W 

-0.  0287 

-0.  0287 

1 

-0.  0287 

0° 

R# 

-0.  239 

-0.  239 

-0.  239 

Mi 

0.  080 

0,  080 

0.  080 

w 

-0.  0149 

-0, 0149 

1 
-0. 0149 

30° 

-0.271 

-0.  271      ; 

-0.271 

•  M6 

0.048   - 

0.  048 

0.  048      j 

w 

0. 0143 

0. 0143 

0.0143 

60° 

"♦ 

-0.  358 

.  -0. 358 

-0.  358 

fi* 

-0.  040 

-0.  040 

-0.  040 

w 

0.  0298 

0.  0298 

-0.  0298 

o 
90 

-0. 413 

-6.413 

-0.413 

M6 

-0.  095 

-0.  095 

-0.  095 
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CHAPTEK  VII 
COMPA  BISON  AND  DISCUSSION  OF  THE  RESULTS 

In  order  to  provide  a  standard  to  which  the  results  of  the 
Fourier  series  and  finite  element  methods  can  be  compared,    the 
ring  problem  is  solved  analytically  for  N.    and  Mi   in  Appendix  A 
for  the  three  load  conditions.      The   results  of  this  analytical 
solution  are  non-dimensionalized  and  are  entered  in  Tables  XII, 
XIII  and  XIV  with  the  results  from  the  Fourier  series  and  finite 
element  methods.     The  Fourier  series  values  in  these  three 
tables  are  the  converged  values  of  w,    N,  and  Mi,   while  the  finite 
element  values  are  the  same  quantities  obtained  by  using  fi  —  7.  5 
degrees.     A   comparison  of  the  results  entered  in  Tables  XII, 
XIII  and  XIV  shows  essentially  zero  error.     This  is  well  within 
normal  engineering  tolerances  for  structural  analysis  and  indicates 
that  under  proper  conditions  either  of  the  two  methods  is  capable 
of  satisfactory  accuracy. 

Since  accuracy  alone  does  not  provide  a  sufficient  basis  for 

choosing  one  method  over  the  other,    the  convenience  or  ease  of 

usage  of  each  method  must  be  examined.      Ease  of  usage  must  be 

considered  within  the  framework  of  the  load  applied,    the  accuracy 

desired,    the  equipment  available,    and  the  assumptions  made.     A 

comparison  of  the  ease  of  usage  of  the  two  methods  within  the 

framework   of     the  load  applied  reveals  that  the  Fourier  series 

method  has   distinct  advantages  over  the  finite  element  method  when 

solving  the  ring  problem  for  a  load  that  is  applied  over  the  surface 

of  the  ring  in  such  a  manner  that  P_  is  a  function  that  is  easily 
6  m  } 

evaluated  for  changing  m  and  (j).     There  are  two  reasons  for  this. 
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TABLE  XII 
RESULTS  FOE  CASE   1  LOAD 


1  ♦ 

Fourier  Series 

Finite  Element 

A  nalytical 

w 

-0.  0744 

-0. 0744 

- 

0° 

N 

0.  000 

0.  000 

0.  000 

0.  318 

0.  318 

0.318 

\v 

-0. 0334 

-0.  0334 

- 

30° 

^ 

-0.  250 

-0.  250 

-0.250       • 

M* 

0.  068 

0.  068 

0.  068 

w 

0.  0364 

0.  0364 

- 

60° 

"♦ 

-0.433 

-0.433 

-0.  433 

-0.  115 

-0.  115 

-0.  115 

w 

0.0683 

0.  0683 

- 

90° 

"♦ 

-0.  500 

-0. 500 

-0.500 

"♦ 

-0.  182 

-0.  182 

-0.  182 

49 


TABLE  XIII 
RESULTS  FOB  CASE  2  LOAD 


* 

Fourier  Series 

Finite  Element 

Analytical 

w 

-0.  0596 

-0. 0596 

- 

0° 

s» 

-0. 128 

-0. 128 

-0.  128 

4> 

0.  190 

0.  190 

0.  190 

w 

-0.  0287 

-0.  0287 

- 

30° 

N«> 

-0.  239 

-0.239 

-0.  239 

M^ 

0.  080 

0.  080 

0.  080 

w 

0.  0298 

0.  0298 

- 

60° 

»t 

-0.413 

-0.413 

-0.413 

M4 

-0.  09  5 

-0.  095 

-0.  095 

w  ,  " 

0.  0574 

0. 0574 

- 

90° 

% 

-0.477 

-0.  477 

-0. 477 

Mf 

-0.  159 

-0. 159 

-0.  159 
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TABLE  XIV 
RESULTS  FOR  CASE  3  LOAD 


4> 

Fourier  Series 

Finite  Element 

A  nalytical 

w 

-0. 0287 

-0. 0287 

"' 

0° 

^ 

-0. 239 

-0. 239 

-0. 239 

"♦ 

0.  080 

0.  080 

-        '   0f-080 

w 

-0.  0149 

-0.0149 

- 

30° 

s 

-0.  271 

-0. 271 

-0.  271   ' 

M^ 

0.  048 

0.  048 

0.  048 

w 

0.  0143 

0. 0143 

60° 

st 

-0.  358 

-0.  358 

-0.358 

u^ 

-0.  040 

-0.  040 

-0.  040 

W 

0.  0298 

0.  0298 

- 

90°' 

4> 

-0.413 

-0.413 

-0.  413 

M^ 

-0.  095 

-0.  095 

-0.  095 
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The  first  is, -that  the  series  is  easily  evaluated  when   P       is  a 
simple  function.     The  second  is  that  for  distributed  loads  the 
finite  element  method  required  initial  fixed  end  conditions  for 
those  elements  containing  the  loads.     Unless  these  initial  end 
conditions  are  available  from  a  handbook  or  other  source.they 
must  be  calculated  using  energy  principles  or  some  other  method. 
As  can  be  seen  in  Appendix  B,    calculation  of  these  initial  end 
conditions  may  be  laborious  and  time  consuming  and  may  re- 
quire a  high  degree  of  initial  accuracy  for  the  results  to  be 
useful.      Even  if  the  equations  for  calculating  the  initial  end 
conditions  are  available  in  a  handbook,    they  may  be  lengthy 
and  may  also  require  a  high  degree  of  initial  accuracy. 

A   situation  in  which  the  finite  element  method  is  superior, 
in  the  framework  of  the  load  applied,    is  one  where  the  load  con- 
ditions lead  to  either    \X0|  —  0  initial  end  conditions  or  fixed  end 
conditions  that  are  available  and  easily  calculated.      This  is 
especially  significant  if  several  such  loading  conditions  are 
to  be  considered.     All  that  is  required  for  a  solution  in  such 
a  situation,    once  a  computer  program  has  been  written,    is  to 
change  the  program  inputs  to  reflect  the  changes  in[x}    , 
\X  \   and    Jul  caused  by  the  load  change  and  to  run  the  program 
for  each  loading  condition.      If  the  Fourier  series  method  is 
used  in  this   situation,    each  load  change  requires  an  integration 
to  obtain  P_j    and  a  different  series  must  be  calculated  for  each 
loading  condition. 

As  pointed  out  earlier  in  this  chapter,    both  methods  are 
capable  of  giving  excellent  accuracy  under  certain  conditions. 
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For  the  Fourier  series  method,    accura>  \  -*d  only  b}    the 

number  of  terms   included  :r.  the  summation,  pre-  ess.      For  dis- 
placemer  t s  .    Table  '1  and  IV  show  that  very  few  terms  are 

required  for  a  high  degrue  of  accuracy  especially  when  a  distri- 
buted load  Is  applied.      Using  only  4  terms   in  any  of  the  summations 
for  w,    including  that  for  the  point  load,    gives  a  maximum  error 
of  0.  02  per  cent.      For  N^  and  Ma  the  situation  is   not  quite  as 
favorable,    but  at  the  worst.    4  terms   give  a  maximum  error  of 
3.  5  per  cent.     If  this   degree  of  accuracy  is   satisfactory,    the 
Fourier   series  method  must  be  considered  easier  to  use  than  the 
element  method. 
For  the  finite  element  method,    all  element  sizes  considered 
gave  essentially  the  same  results  for  all  three  load  conditions. 
However,    for  load  conditions  o'her  than  the  ones  i     ••ordered  here, 
the  aua'vst  may  be  required  to  use  elements  with  (i  less  than  7.  5 
degrees.      There  appear  to  be  two  limiting  factors  <  r   the  minimum 
siz<=   of  elements  that    can  be  used.    The  first  limiting  factor  is 
puter  capacity.     As   the  el  ts  are  mad-     -.-nailer,    more 

nodes  are  created,    and  hence  larger  matrices  are  required.      When 
dealing  with  full  scale  structures,    the  limit  of  computer  capacity 
can  be  reached  prior  to  reaching  the  optimum  element  size  for 
accuracy.      The  oth'.r  limiting  factor  is  related  to  computer 
capacity  bit  a?  rom  a  charact<  curved  element 

ffness  matrix   rather  than  from  the  number  of  elements  used. 
The  characteristic  referred  to  is  the  tendency  for  some  of  the 
numb'  matrices,    particularly  on  the  main  diago- 

nal ,    to  have  a   greater  dispax  size  as  (I    Le  <l<  creased.     Tables 

I  and  VII  show  this  to  a  marked  degree.      For  instance,    the 
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ratio  of  the  stiffness  coefficient  (6,  6)  to  coefficient  (1,  1)  in 

K,,  '•  for    (3   —  30  degrees   is   17/17,  243,    whereas  the  ratio  of 
the  same  two  elements  is  69/  18,  665,  685  for  ft  =   7.  5  degrees. 
Since  the  computer  retains  only  a  certain  number  of  significant 
figures  while  manipulating  these  matrices,    a  point  will  be 
reached  where  round-off  error  in  computer  operations  will 
negate  any  gain  in  accuracy  obtained  by  using  smaller  elements. 
This  problem  is  alleviated  to  a  degree  by  using  double-precision 
techniques,   as  used  in  this  thesis;   but  again  since  double- 
precision  operations  require  more  computer  storage  space,    the 
computer  capacity  itself  places  a  limit  on  this  approach. 

The  equipment  available  to  the  analyst  has  a  significant 
effect  on  the  ease  of  usage  of  the  two  methods.     Indeed,    in  the 
case  of  the  finite  element  method,    lack  of  a  digital  computer 
almost  prohibits  use  of  the  method  on  the  basis  of  effort  required. 
Even  if  a  digital  computer  is  available,    a  desk  calculator  and  a 
set  of  mathematical  or  trigonometric  tables  with  up  to  ten  signifi- 
cant figures  is  desirable  if  the  nature  of  the  load  is  such  that 
fixed  end  conditions  must  be  determined.      On  the  other  hand, 
while  a  desk  calculator,    and  even  a  digital  computer,    may  be  used 
to  advantage  in  the  Fourier  series  method,    neither  are  required.-' 
A   slide-rule,    a  table  of  integrals,    and  normal  trigonometric 
tables  are  sufficient  equipment  for  application  of  the  Fourier 
series  method. 

Finally,    the  ease  of  usage  must  be  considered  within  the 
framework  of  the  theory  employed.      The  assumptions  made  when 
applying  the  two  methods  of  solution  to  the  ring  can  have  an  effect 
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on  both  the  accuracy  obtained  and  the  e,   5      . M   u  s  igc. 
governing  differential  equations  are  derived,    and  the  Fourier 
series   solution  developed  through  Equation  (4-6a),    (4-6b)  and 
(4-6c)  without  making  any  assumptions   regarding  the  extension 
of  the  middle  surface  of  the  ring.     If  the  assumption  is  made 
during  the  derivation  of  these  differential  equations  that  the 
middle  surface  is  inextensional,    that  is  £  ^-  0    when  z  -  0,    then 
-r-r —  z  -w  and  the  m  ~    0  ,    or  extensional,    terms  do  not  appear  in 
the  expressions  for  w  and  M.  .  *      As   seen  in  Table  I,    no  large 
error  is  introduced  in  w  and  M.   by  this  assumption.      However, 
if  this  inextensional  form  of  w  is  used  in  Equation  (2-6a)  to 
obtain  N.,    the  resulting  expression  is  missing  the  leading,    or 
m   —    0,    term,    and  the  values  for  N.   calculated  from  this  shortened 
expression  are  in  unacceptable  error.      This  problem  may  be 
circumvented  by  using  the  equilibrium   equation,    Equation  (2-2a), 

to  calculate  N,.  The  procedure  used  is  to  substitute  the  inex- 

f 

tensional  form  of  M^,    obtained  from  Equation  (2-6b),    into 
Equation  (2-2a)and  integrate.      The  constant  of  integration  is 
obtained  from  equilibrium  considerations  at  one  of  the  planes  of 
symmetry.      The  resulting  expression  for  N,    is  identical  to 
Equation  (4-4b)  and  thus   gives  accurate  results.      On  the  other 
hand,    in  the  finite  element  method,    an  error  will  occur  in  the 


:  The  governing  equations  for  the  modes  m^2  identically 
satisfy  the  inextensional  assumption.      Thus,    these  nodes  are 
not  affected  by  the  assumption. 

This  is  analgous  to  determining  Qa  from  the  equilibrium 
equation,    Equation  (1-la). 
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solution  for  w  if  the  inextensional  assumptioa  is  m-de  ar.d 
distributed  loads  are  applied  over  almost  all  the  ring  surface. 
The  reason  for  this  is  that  a   distributed  load  acting  on  the  surface 
of  a  curved  element  has  a  greater  tendency  to  compress  or  ex- 
tend the  middle  surface  than  does  a  concentrated  load.      Errors 
are  caused  if  the  energy  of  this  middle  surface  compression  or 
extension  is  neglected,   as  it  is  in  the  derivation  of  the  curved 
element  stiffness  matrix.      Unlike  the  Fourier   series  method, 
there  is  no  means  of  introducing  this  extensional  behavior  into 
the  finite  element  method  once  the  stiffness  matrix  is  derived. 
The  only  recourse  is  to  select  a  t/a  ratio  and  re-derive  the 
curved  element  stiffness  matrix  including  the  axial  strain  energy 
term  in  the  calculations. 
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CHAPTER  VIII 
CONCLUSIONS  AND  RECOMMENDATIONS 

The  primary  advantage  of  the  Fourier  series  method  for 
analyzing  ring  structures  is  that  it  gives   reasonably  accurate 
answers  with  a  minimum  amount  of  effort  using  a  slide  rule 
and  normally  available  integral  and  trigonometric  tables. 

The  primary  advantage  of  the  finite  element  method  lies 
in  its  ability  to  accurately  solve  several  different  loading  con- 
ditions with  a  minimum  amount  of  effort.     The  finite  element 
method  has  the  disadvantage  of  requiring  the  calculation  of 
fixed  end  conditions  for  certain  types   of  loads. 

A   recommendation  for  further  work  is  the  derivation  of  a 
stiffness  matrix  for  a  curved  element  that  includes  the  effects 
of  axial  strain  energy.     Also  recommended  is  a  matrix  error 
analysis  on  both  the  extensional  and  inextensionai  forms  of  the 
curved  element  stiffness  matrix  for  various  values  of  the 
central  angle. 
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APPENDIX  A 
ANALYTICAL  SOLUTION 

The  ring  problem  is  solved  analytically  for  Ni  and  Mx  for 
all  three  load  conditions  to  obtain  a  standard  to  which  the  Fourier 
series  and  finite  element  results  can  be  compared.     The  ana- 
lytical solution  is  based  on  the  symmetric  character  of  the  loads 
applied,    the  principle  of  superposition  and  Castigliano's  Theorem. 
The  use  of  superposition  in  the  solution  of  thin  shell  problems  is 
discussed  in  detail  in  Reference  5.     Figure  A    1  (I)  shows  the 
entire  ring  with  an  applied  load  of  P  pounds  uniformly  distri- 
buted over  each  of  the  two  surface  areas  subtended  by  the  angle 
2©<  .     Figures  A    1  (II),    (III)  and  (IV)  show  the  portions  of  the 
ring  that  are  superimposed  to  obtain  a  solution. 

Considering  the  quadrant  of  the  ring  for   0  —  #vf  tT/ 2 ,      the 
boundary  conditions  plus  equilibrium  for  Figure  A    1  (II)  require 
that 

.     (N^)n    =  -    P/i~  (A-la) 

and 

(^)n  =    -    0  (A -lb) 

for0^4)-°<.     Similarly 

<VlII    =       H2COS     ♦  (A-Za) 

and 

(Mi)         =       M2  +    aH2    (cos     (J)     -    cosx)  (A -2b) 

for  Figure  A    1  (III). 
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FIGURE  A    1 
RING  SEGMENTS  USED  IN  ANALYTICAL  SOLUTION 
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In  Figure  A    1  (IV) 

Psin  o<  sin     (|) 

>  'IV 

and 


(NJn/  2o<  ^  -3a) 


a  Psinoi  ( 1  -sin 
(VlV    =       Mi+  2^ 


are  written  for  o<   ^  (j)  ^   TT/2.     Since  equilibrium   gives 

aP  (l-sino<  ) 


M,  —      M2  -  aH£    cos  <=< 


M 


,   can  be  substituted  in  the  expression  for  (M(k)TV  to  obtain 

i\a    \                   \a              u                 ^1    aP  (l-sino<sin  (D)  /A     ou » 
(M<(>TV      ~      M2  -  a  H2  cos  c<  -|-  2o< "  (A  -3b) 


When  $ 

=r  <=*; 

\ 

Mz 

so 

that 

Equation 

(A- 

■3b) 

yields 

H2 

= 

Pcosot 
2  o< 

This  expression  is  substituted  into  Equations  (A -2a)  and  (A -2b), 
and  the  resulting  equations  are  superimposed  on  Equations  (A -la) 

and  (A  -lb)  to  obtain 

p 

Ni,  -         r-^    (coscx  cos      <|)-1)  (A-4a) 


and 

<P  2  2  <* 


M,  =  M    +     - — ?  (cos   <()  -  cos°<)  (A-4b) 


for  0  -  ty-<*. 

Substituting  the  same  expression  into  Equation  (A  -3b)  yields 

M.  =         M.    +     a  F  sin    (sin-<-     sin  0)-  (A -5) 

9  2  ^ 

where  c*  *    (p    *      TT/2. 

Castigliano's  Theorem  is  applied  in  the  form   of 
L  ^ty/  =0  to  obtain  the  internal  moment  ML.      The  total 

strain  energy,    U  ,    of  the  three  segments  is  written  as 

Ut  =        UIIar  UHb '  UIIIa •'■  UIIIb+uIVa  fUrVb 
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where  subscripts  a  and  b  stand  for  axial  and  bending  strain 
energy  respectively.     The  terms  Utt    ,    Utt,  ,    U-q^,    and  Ujya  need 
not  be  considered  when  applying  Castigliano's  Theorem  to  obtain 
M?  since  these  terms  are  either  zero  or  are  not  functions  of  Mi. 


Therefore 


^u 


y^Mi 


=   0 


becomes 


M 


2.  EI 


^-ad4> 


+■ 


-1T/z     .1 


<-  •'O 


=    0 


(A-6) 


where  the  expressions  for  M     are  obtained  from  Equations  (A -4b) 
for  the  first  integral   and  (A -5)  for  the  second  integral.  Performing 
the  partial  differentiation  and  integration  with  the  appropriate 
limits  as  indicated  in  Equation  (A  -6)  gives 


M,  =  - 


1T 


1T     MY\     <=< 


2.  << 


-J 


(A -7) 


Using  this  expression  for  M     in  Equations  (A -4b)  and  (A -5),    the 
complete  set  of  equations  for  Ni  and  Mi  is  written. 


No=  T^ 


COW    COS  ([)    -J 
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M*  =  aP 


J. 

T 


2°< 


C  Ob^    COS    if       ] 


for   0   ^    f  ^  °<  ,    and 


»  L    -  - 


N<,= 


P  Sm°^   sm  $ 


Z<* 


M*=aP 


r 

a- 

ir 


S  l  vi.  oc     s  \  n  c|- 


^ 


for<X   Zl   <f  —  TT/2.  •  These  equations  are  non-dimension- 

alized  and  evaluated  for  ((1=0,    30,    60  and  90  degrees  for  all 
three  case  loads.     The  results  are  entered  in  Table  AI. 
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TABLE  A  I 
ANALYTICAL  SOLUTION  EESULTS 


<p 

Case   1  Load 

Case  2  Load 

Case  3  Load 

0.  0 

-0.  128   ' 

-0.  238 

0° 

^ 

0.  318 

0.  190 

0.  080 



** 

-0.  250 

-0. 239 

-0.  270 

3-s°; 

M* 

-0.  068 

0.  080 

0.  048 

N* 

-0.433 

-0.413 

-0.358 

60° 

. 

M^ 

-0.  115 

-0.  095 

-0.  040 

N^ 

-0.  500 

-0.477 

-0.413 

90° 

M* 

-0.  182 

-0. 159 

-0.  095 

'- 
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APPENDIX  B 
SOLUTION  FOR   INITIAL  FORCES 

The  techniques  used  to  obtain  the  initial  fixed  end  moments 

and  forces   required  for  the  finite  element  solutions  to  Case  2  and 

Case  3  loads  involve  superposition  and  Castigliano's  Theorem. 

Figure  B   1   shows  the  finite  element,    the  three  loading  conditions, 

the  boundary  forces  and  moments,    and  the  notation  used.      The 

superposition  of  the  forces  and  moments  in  Figure  B   1  (II)  and 

B   1  (III)  is  equivalent  to  the  forces  and  moments  in  Figure  B   1  (I). 

The  sum  of  the  rotations   of  the  ends  of  the  element  in  Figure 

B   1  (II)  and  (III)  must  be  zero.      Thus,    en  -r    9m    =    0. 

However,    since  6      is  zero  under  the  end  conditions  and  loading 
II 

of  Figure  B   1  (II),    0ttj  =   0.      The  analogous   requirement  for  the 
horizontal  displacements  of  the  ends  of  the  element  of  Figure 
B   1  (II)  and  (III)  leads  to 


^        =    ^  (B-l) 


where     ^>_      is  the  horizontal  component  of  the  known  radial  dis- 
placement of  (II),    and     o>rn  ^s  th.e  horizontal  displacement  obtained 
from  an  expression  for  the  strain  energy  of  the  element  in  Figure 
B   1  (III)  using  Castigliano's  Theorem.         ^-p,    is  positive  in  the 
direction  towards  the  centerline,    and        o-nr     *s  Posltive  in  the 
direction  away  from  the  centerline. 
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Letting  o<  be  the  half  angle  of  the  element  and  (j)  be  the 
angle  at  any  point  along  the  element,    the  moment  at  (j)  is 

I  r  oc    rh     -    C  0S«=^  j 


M$  =  M   -aHm'cos  0  -  c 


And  the  axial  force  at  Q  is 


\p  =    Hut  cos  $ 


The  strain  energy,    U    from  both  axial  extension  and  bending 
of  the  total  element  is 


Ut=  2. 


ZEl  I       2AE 


r 

/ 

/ 
'  o 


(B-2) 


Using  Equation  (B-2)  in  Castigliano's  Theorem  gives 


Li 


>>v* 


Substituting  Equation  (B-2)  into  Equation  (B-3),    differentiating, 
integrating  and  applying  limits  yield 


(B-3) 


^m 


a 


^.HnJr 


$m©<-<*  co 


c  cv 


3  vn  *■«* 


+■  «*  cos  <=* 


an 


nr 


M: 


ex 


(B-4) 


+ 

f       J 


6? 


The  horizontal  displacement   o  TT    is 


W= 


a  P  sm 


ex. 


z^V>-HE 


where  f  is  given  by 


(B-5) 


CASE  LOAD 

2 

3 

<3  =      7.5° 

f 

=     8 

f   =    16 

0=  15.0° 

f 

=     4 

f    =      8 

(?  =  30.  0° 

f 

=     2 

f    =      4 
> 

The  factor  f  is  required  in  order  to  give  the  proper  portion  of  the 
total  load  P  over  each  of  the  different  segment  sizes.     Substi- 
tuting Equations  (B-4)  and  (B-5)  into  (B-l)  leads  to 


flH™ 


S  \V\<=<-<X  C 


0S~J    +  — - 


3S\Y\2<*  '     2 


♦£"* 


—  +-    

2-  if 


a  P  sin** 


The  above  equation  can  be  simplified  to  give 

H(TllHJ-^tsvn^-*coH>T- 


3  s  \  Y\  1^ 


+  c^    coS=< 


-hH 


HI 


SI 


Y\  2oc]  P  S  i  h 


oC 


*f 


2.c*f 


(B-6) 


Thus,    H„.can  be  determined  from  Equation  (B-6)  for  a  given  case 
load  and  a  specified  central  angle  for  the  element,  /S  ,    once  t/a  is 
chosen.      The  requirement  to  select  a  specific  t/a  ratio  in  order  to 
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solve  Equation  (B-6)  means  that  the  axial  strain  energy  is  in- 
cluded in  the  expres-sion.      Failure  to  include  the  axial  strain 
energy  by  making  the  assumption  made  in  Reference  2  that  it  can 
be  neglected  for  curved  elements  with  small  t/a  ratios  leads  to 
gross  errors.     This  is  because  the  assumption  made  in  Refer- 
ence 2  is  not  valid  for  curved  elements  with  distributed  loads 
over  large  portions   of  the  element. 

The  equation  for  the  total  normal  forces  N°  acting  on  the 
ends  of  the  element  shown  in  Figure  B  1  (I)  is 


or 


N    =   Nm-  Hmcoso< 


N°=    —  -  Hm  cos~<  (b-7) 


since 

nt  of  1 
force,    Q°,    at  the  ends  of  the  ring  element.     Thus 


The  component  of  H       in  the  radial  direction  is  the  shear 


Q°=  Hm    SIM  (B-8) 

In  order  to  obtain  M    ,    Castigliano's  Theorem,    in  con- 
junction with  the  requirement  9ttt  =   0,  is  used  to  write 

where  U    is  given  by  Equation  (B-2). 
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Differentiation  and  integration  of  this  equation  give 


m0  =  «v(  J^fL -C0H  (B-9j 


Experience  with  Equations  (B-6)  through  (B-9)  shows  that  the 
axial  strain  energy  must  be  included  if  accurate  answers  are  de- 
sired for  N°,    Q°  and  M    .     Thus,    t/a    =    1/10  is  chosen,and 
Equations  (B-6)  through  (B-9)  are  used  to  prepare  the  entires 
in  Table  B  I  for  N°,   and  M     for  Case  2  and  3  loads  and    p     =      30, 
15  and  7.  5  degrees.     In  evaluating  these  equations,    trigonometric 
functions  accurate  to  5  decimal  places  do  not  give  accurate 
answers  for  N    ,  Q     and  M     for  the  values  of    /->    used  in  this 
thesis.     Consequently  the  entries  in  Table  B  I  are  obtained  by 
using  trigonometric  functions  of  10  decimal  place  accuracy. 
Solving  equations  such  as  (B-6)  using   10  decimal  places  is  ob- 
viously beyond  the  capacity  of  either  slide  rule  or  practical  hand 
calculations.     Even  with  an  electronic  desk  calculator,    solving 
such  equations  is  laborious  and  time  consuming. 
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N 


M 


TABLE  B  I 
FIXED  END  RESULTS 
CASE  2  LOAD 


7.5" 


0.  00181929P 
0.  06247015P 
0.  00136328Pa 


15° 
0.  01286764P 
0.  12402482P 
0.  00541780Pa 


30° 


0.  12707678P 
0.  22182251P 
0.  01944670Pa 


CASE  3  LOAD 


7.5° 

15° 

3  0° 

N° 

0.  00090964P 

0.  00643382P 

0.  06353839P 

Q° 

0.  03123508P 

0.  06201237P 

0.  1109H25P 

M° 

0.  00068164Pa 

0".  00270890Pa 

0.  00972335Pa 
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APPENDIX  C 

FORTRAN  IV  PROGRAMS 

TABLE  C  I 

DESCRIPTION  OF  SUBROUTINE  "DSIMQ1 


PURPOSE 

OBTAIN  SOLUTION  OF  A  SET  OF  SIMULTANEOUS  LINEAR 
EQUATIONS,   AX=B 

USAGE 

CALL  "DSIMQ"  (A,B,N,KS) 

DESCRIPTION  OF  PARAMETERS 
A  AND  B  MUST  BE  REAL*8 
A    -  MATRIX  OF  COEFFICIENTS  STORED  COLUMNWISE. 

THESE  ARE  DESTROYED  IN  THE  COMPUTATION. 

THE  SIZE  OF  MATRIX  A  IS  N  BY  N. 
B  -    VECTORS  OF  OR IGINA  L  CONSTANTS  (LENGTH  N). 

THESE  ARE  REPLACED  BY  FINAL  SOLUTION  VALUES, 

VECTOR   X. 
N  -    NUMBER   OF  EQUATIONS  AND  VARIABLES 
KS  -OUTPUT  DIGIT 

0  FOR  A  NORMAL  SOLUTION 

1  FOR  A  SINGULAR  SET  OF  EQUATIONS 

REMARKS 

MATRIX  A   MUST  BE  GENERAL. 

IF  MATRIX  IS  SINGULAR,    SOLUTION  VALUES  ARE  MEAN- 
INGLESS.    AN  ALTERNATIVE  SOLUTION  MAY  BE  OB- 
TAINED BY  USING  MATRIX  INVERSION  (MINV)  AND 
MATRIX  PRODUCT  (GMPRD). 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 

NONE 

METHOD 

METHOD  OF  SOLUTION  IS  BY  ELIMINATION  USING 
LARGEST  PIVOTAL  DIVISOR.     EACH  STAGE  OF  ELIMINA- 
TION CONSISTS  OF  INTERCHANGING  ROWS  WHEN  NECES- 
SARY TO  AVOID  DIVISION  BY  ZERO  OR  SMALL  ELEMENTS. 
THE  FORWARD  SOLUTION  TO  OBTAIN  VARIABLE  N  IS 
DONE  IN  N  STAGES.     THE  BACK  SOLUTION  FOR   THE 
OTHER  VARIABLES  IS  CALCULATED  BY  SUCCESSIVE 
SUBSTITUTIONS.     FINAL  SOLUTION  VALUES  ARE 
DEVELOPED  IN  VECTOR   B,    WITH  VARIABLE   1  IN  B  ( 1 ), 

VARIABLE  2.  IN  B(2) ,   VARIABLE  N  IN  B(N). 

IF  NO  PIVOT  CAN  BE  FOUND  EXCEEDING  A   TOLERANCE 
OF  0.  0,    THE  MATRIX  IS  CONSIDERED  SINGULAR  AND 
KS  IS  SET  TO  1.     THIS  TOLERANCE  CAN  BE  MODIFIED 
BY  REPLACING  THE  FIRST  STATEMENT. 
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TABLE  C  II 
LIST  OF  SYMBOLS  USED  IN  SUBROUTINE  "RELM" 


Computer 
Coded  Name 

Definition 

DBETA 

/3  ,    element  central  angle. 

DESM 

Element  stiffness  matrix      K  «  j  . 

DESM    (I, 

J) 

Stiffness  influence  coefficient  in 

DLA 

a            (corresponding  notation  in 

DLB 

b             Reference  3). 

DLC 

c 

DLD 

d 

DLE 

e 

DBA 

A 

DBB 

B 

DBC 

C 

DBD 

D 

DBE 

E 

DBF 

F 

DBG 

G 
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TABLE    CIII 
SUBROUTINE    "RELM" 

SUBROUTINE    RELM(  DBETA ,DESM, DL A, DLB, DLC, OLD, DLE , DBA, 
1DBB,DBC,DBD.DBE,DBF,D8G) 

DOUBLE    PRECISION    DLA, DLB, DLC, DLD, DLE, DBA,DBB,DBC,DBD, 
1DBE,DBF,DBG 

DOUBLE    PRECISION    DBETA 

REAL*8    DESM(6,6) 

DLA=OBETA-DSIN< DBETA) 

DLB=DCOS< DBETA) ♦(< OS IN< DBETA)  )**2 ) /2.0D0-1 .ODO 

DLC=   (3.0D0*DBETA)/2.0D0-2.0D0*DSIN(DBETA)+ 
l(DSIN(2.0DO*DBETA) )/4.0DO 

DLC=DBETA/2.0D0-(DSIN(2.0D0*DBETA) )/4.0DO 

DLE=DCOS( DBETA )-l. ODO 

DBA=(DLE**2)/DBETA-DLD 

DBB=OLB-(DLA*DLE)/ DBETA 

D8C=0LA*DLD-DLB*DLE 

DBD=(DLA**2) /DBETA-DLC 

D3E=DLC*DLE-DLA*DLB 

D3F=(DLB**2)-(DLC*DLD) 

DBG=DLB*(DLB-(2.0DO*DLA*DLE)/DBETA)+DLC* 
1 ( -OLD* (( DL E**2 I / DBETA )  )♦< ( DLA**2 )*DLD) /DBETA 

DESM(1,1)=DBA/DBG 

DESMU,2)=DBB/DBG 

DESM(1,3)=DBC/(08ETA*DBG) 

DESM(2,1)=DBB/DBG 

DESM(2,2)=DBD/DBG 

DESM(2,3)=DBE/(DBETA*DBG) 

DESM<3,1)=DBC/(DBETA*DBG) 

DESM(3,2)=DBE/(DBETA*DBG) 

DESM(3,3)=DBF/(DBETA*DBG) 

DESM(4,1)=-(DBA*DC0S(DBETA)+DBB*DSIN(DBETA) I /DBG 

DESM<4,2)=-    (DBB*DCOS(DBETAKDBD*DSIN(  DBETA)) /DBG 

DESM(4,3)=-       (DBC*OCOS(DBETA)+DBE*DSIN(DBETA) )/ 
l(DBETA*DBG) 

DESM(5,1)=  <DBA*DSIN(D8ETA)-DBB*DCOS( DBETA) )/DBG 

DESM(5,2)=       <DBB*DSIN(DBETA)-DBD*DCOS( DBETA)) /DBG 

DE SMC  5,3)=       <DBC*DSIN(DBETA)-DBE*DCOS(DBETA) )/ 
1(DBETA*DBG) 

DESM(6,l)=  (DBA*(DC0S<D8ETA)-1.0D0H-0BB* 

KDSIN(DBETA)  )-DBC/DBETA)  /DBG 

DESM(6,2)=       (DBB*(DCOS(DBETA)-1.0DO)+DBD* 
IDS IN( DBETA J-DBE/DBETA) / DBG 

DESH(6.3)=       (DBC*(DCOS(DBETA)-1.0DO)+DBE* 
1DSIN(DBETA)-DBF)/(DBETA*DBG) 

OESM(4,4)=DESM(l.l) 

DESM<4,5)=-DESM(1,2) 

DESM(4,6)=DESM(1,3I 

DESM(5,4)=-0ESM(2,1) 

DESM(5,5)=DESM<2.2) 

DESM(5,6)=-DESM(2,3) 

DESM(6,4)=DESM(3,1) 

DESM(6,5)=-DESM(3,2) 

DESM(6,6)=DESM(3,3) 

DESM(l,4)=OESM(4,l) 

OESM(l,5)=DESM(5,l) 

DESM(1,6)=DESM(6,1) 

DESM(2,4)=DESM(4,2) 

DESM<2,5)=DESM(5,2) 

DESM(2,6)=DESM(6,2) 

DESM<3,4)=DESH<4,3) 

DESM<3.5)=DESM(5,3) 

DESM(3,6)=DESM(6,3) 

RETURN 

END 
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TABLE  C  IV 
LIST  OF  SYMBOLS  USED  IN 
PROGRAM  "BETA   30" 


Computer 
Coded  Name 

A 

A  (I,  J) 

B  (I) 

B  (I) 
DBETA 
DESM 
DESM  (I,  J) 

PI 

DM  (I) 
DN  (I) 


Definition 

Reduced  stiffness  matrix     Kp   . 

Stiffness  influence  coefficient 
in    J  K/5 '  . 

Prior  to  calling  "DSIMQ"  B(I)  = 

\  M     -      fx°|]  • 

After  calling  "DSIMQ"  B    (I)  =   [uf\ 
/3     ,    element  central  angle. 
Element  stiffness  matrix      |K« 

Stiffness  influence  coefficient 

in    !K, 


r\ 


tr 


Moment  at  node  (I). 
Normal  force  at   node  (I). 
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TABLE    C5T 

PROGRAM    "BETA30" 

PROGRAM  BETA  30 

DOUBLE  PRECISION  DBETA. PI , DLA, DLB, DLC, DLD, OLE , DBA, 
lDBB,DBC,DBD,DBEf DBF, DBG 

91  FORMAT  ( 1H0,1P6D20.10) 

92  FORMAT  ( 1H0 , 1P5D19 . 10  ) 

93  FORMAT  ( 1H0 , 1P7D18 . 10 ) 

94  FORMAT  ( 1H0 , 1PD20. 10 ) 

11  FORMAT  (1H0,1P5D20.10) 

12  FORMAT  (1H0, 1PD20.10) 

13  FORMAT  <1H0,1P6D20.10) 

14  FORMAT  < 1H0,1P9D15.7) 

15  FORMAT  (1H1, •DISPLACEMENTS1 ) 

16  FORMAT  (1HO,»FORCES  AND  MOMENTS') 

17  FORMAT  (lHO.'CASE  2  BETA=30  DEGREES1) 

22  FORMATdHO,  1PD20.10) 

23  F0RMATU2) 

24  F0RMATUH0.1PD20.10) 

PI =3.1415926535  8979300 

DBETA=PI/6.0D0 

WRITE  (6,94)  DBETA 

REAL*8  DESM(6,6) 

CALL  RELM  ( DBETA, DESM, DLA, DLB, OLC, DLD, OLE, DBA, 
10BB,DBC,DBD,DBE,DBF,DBG) 

WRITE  (6,92)  DLA, DLB. DLC, DLD. OLE 

WRITE  (6,93)DBA,DBB,DBC,DBD,0BE,DBF,DBG 

WRITE  (6,91)  ((DESM( I , J ) , J=l .6 ) , 1= 1, 6) 
F0RBETA=30DEGREES 

REAL*8A<8.8) 

DO  110  1=1,8 

00  111  J=l,8 

A< I, J) =0.000 
111  CONTINUE 
110  CONTINUE 

A(l,l)=DESM(2,2) 

A(2,1)=DESM(4,2) 

A(3,1)=DESM(5,2) 

A<4,1)=DESM(6,2) 

A(i,2)=DESM(2,4) 

A(2,2)=0ESM(4,4)*DESM(1,1) 

A(3,2)=DESM(5,4)*DESM<2,1) 

A(4,2)=DESM(6,4)+DESM(3,1) 

A(5,2)=DESM(4,1) 

A(6,2)=0ESM(5,1) 

A(7,2)=DESM(6,1) 

A(1,3)=DESM(2,5) 

A(2f3)=DESM(4,5)+DESM(l,2) 

A(3,3)=DESM<5,5)+DESM(2,2) 

A(4,3)=DESM(6,5)>DESM(3,2) 

A(5,3)=DESM(4,2) 

A<6,3)=DESM<5,2) 

A(7,3)=DESM(6,2) 

A(1,4)=0ESM(2,6) 

A(2,4)=DESM(4,6)+DESM(1,3) 

A(3,4)=DESM(5,6)+DESM(2,3) 

A<4f4)=DESM(6,6)+DESM(3,3) 

A(5,4)=DESM(4,3) 

A(6,4)=DESM<5,3) 

A(7,4)=DESMC6,3) 
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